﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using T_BPI.Classes.Base;

namespace T_BPI.Classes.Model
{
    public class GMDH
    {
        int m;                 // число признаков
        int N;                 // число образцов выборки (N = NA + NB)
        int NA;
        int NB;
        double percentage;     // процент отбираемых "лучших" моделей
        List<Row> rows = new List<Row>();

        public GMDH(int m, int N, int NA)
        {
            this.m = m;
            this.N = N;
            this.NA = NA;
            this.NB = N - NA;
            percentage = 0.4;
        }

        int comb2(int n)
        {
            if (n == 0 || n == 1) return 0;

            int comb = 0;
            for (int i = 0; i < n; i++)
                for (int j = i + 1; j < n; j++)
                    comb++;
            return comb;
        }

        /* Решение СЛУ методом Гаусса (Жордана) для матрицы [n][n+1]
где, посл столбец - столбец своб членов */
        double[] solveLinearEquasion(double[,] A_matr, double[] b_matr, int n)
        {
            int i, j, k;
            double u;
            double[] x = new double[n];

            // A_matr[n][n]
            // B_matr[n]
            // z[n][n+1] = A_matr B_matr

            double[,] z = new double[n, n + 1];
            for (i = 0; i < n; i++)
            {
                for (j = 0; j < n; j++)
                {
                    z[i, j] = A_matr[i, j];
                }
            }
            for (i = 0; i < n; i++)
            {
                z[i, n] = b_matr[i];
            }

            // приведение матрицы к н. треуг. виду (по Гауссу)
            for (j = 0; j < n + 1; j++)
                for (i = j + 1; i < n; i++)
                {
                    u = -z[i, j] / z[j, j];
                    for (k = j; k < n + 1; k++)
                        z[i, k] = z[i, k] + u * z[j, k];
                }

            // приведение правой матрицы к в. треуг. виду (по Жордану)
            for (j = n - 1; j > 0; j--)
                for (i = j - 1; i >= 0; i--)
                {
                    u = -z[i, j] / z[j, j];
                    for (k = j; k < n + 1; k++)
                        z[i, k] = z[i, k] + u * z[j, k];
                }

            for (i = 0; i < n; i++)
            {
                x[i] = z[i, n] / z[i, i];
            }

            return x; // вектор решений СЛУ
        }

        public void training(double[,] input_data, double[] y0)
        {
            // формирование рядов
            int F = m;
            double prev_error = Double.MaxValue;

            double[,] inputs = new double[N, m];


            for (int i = 0; i < N; i++)
            {
                for (int j = 0; j < m; j++)
                {
                    inputs[i, j] = input_data[i, j];
                }
            }
            double[,] outputs;
            int row_number = 0;

            while (row_number < 50)
            {
                int l = 0;
                row_number += 1;

                int L = comb2(F);    // число моделей текущего ряда = C(F, 2)
                Row row = new Row(L);    // создаем новый ряд

                int[,] combinations = new int[L, 3]; // (model_num, i, j)

                l = 0;
                for (int i = 0; i < F; i++)
                    for (int j = i + 1; j < F; j++)
                    {
                        combinations[l, 0] = l; // model_num
                        combinations[l, 1] = i;
                        combinations[l, 2] = j;
                        l++;
                    }

                // смотрим на обучающую выборку и вычисляем весовые коэффициенты
                for (l = 0; l < L; l++)
                {
                    int model_num = combinations[l, 0];
                    int i = combinations[l, 1];
                    int j = combinations[l, 2];

                    // по всем сочетаниям i и j
                    double[,] A_matr = new double[Constants.WEIGHTS_CONST, Constants.WEIGHTS_CONST];
                    for (int it1 = 0; it1 < Constants.WEIGHTS_CONST; it1++)
                        for (int it2 = 0; it2 < Constants.WEIGHTS_CONST; it2++)
                            A_matr[it1, it2] = 0.0;

                    double[] b_matr = new double[Constants.WEIGHTS_CONST];
                    for (int it = 0; it < Constants.WEIGHTS_CONST; it++)
                        b_matr[it] = 0.0;

                    for (int k = 0; k < NA; k++)
                    {
                        double xi = inputs[k, i];
                        double xj = inputs[k, j];

                        A_matr[0, 0] += 1.0;
                        A_matr[1, 0] += xi;
                        A_matr[2, 0] += xj;
                        A_matr[3, 0] += xi * xj;

                        A_matr[0, 1] += xi;
                        A_matr[1, 1] += xi * xi;
                        A_matr[2, 1] += xi * xj;
                        A_matr[3, 1] += xi * xi * xj;

                        A_matr[0, 2] += xj;
                        A_matr[1, 2] += xi * xj;
                        A_matr[2, 2] += xj * xj;
                        A_matr[3, 2] += xi * xj * xj;

                        A_matr[0, 3] += xi * xj;
                        A_matr[1, 3] += xi * xi * xj;
                        A_matr[2, 3] += xi * xj * xj;
                        A_matr[3, 3] += xi * xi * xj * xj;

                        b_matr[0] += y0[k];
                        b_matr[1] += y0[k] * xi;
                        b_matr[2] += y0[k] * xj;
                        b_matr[3] += y0[k] * xi * xj;
                    }

                    double[] weights = solveLinearEquasion(A_matr, b_matr, Constants.WEIGHTS_CONST);

                    // добавляем модель в ряд
                    row.models[l] = new Model(model_num, i, j);
                    row.models[l].setWeights(weights);
                }

                // вычисление ошибки для каждой модели
                for (l = 0; l < L; l++)
                {
                    int model_num = combinations[l, 0];
                    int i = combinations[l, 1];
                    int j = combinations[l, 2];

                    row.models[model_num].CR = 0.0;
                    // по всем образцам из проверочной выборки
                    for (int k = NA; k < N; k++)
                    {
                        double xi = inputs[k, i];
                        double xj = inputs[k, j];
                        row.models[model_num].CR += (row.models[model_num].f(xi, xj) - y0[k]) * (row.models[model_num].f(xi, xj) - y0[k]);
                    }
                    row.models[model_num].CR /= NB;
                }

                row.countCR();

                bool changes = true;
                while (changes)
                {
                    changes = false;
                    for (l = 0; l < L - 1; l++)
                    {
                        if (row.models[l].CR > row.models[l + 1].CR)
                        {
                            swapModelsInArray(row.models, l, l + 1);
                            changes = true;
                        }
                    }
                }

                F = (int)this.percentage * L; // число отбираемых лучших моделей
                if (F > Constants.MAX_MODELS_NUM) F = Constants.MAX_MODELS_NUM;

                outputs = new double[N, F];
                for (int k = 0; k < N; k++)
                    for (l = 0; l < F; l++)
                    {
                        double xi = inputs[k, row.models[l].input1];
                        double xj = inputs[k, row.models[l].input2];
                        outputs[k, l] = row.models[l].f(xi, xj);
                    }

                if (prev_error < row.CR)
                    return;
                prev_error = row.CR;
                rows.Insert(rows.Count, row);

                if (F <= 1)
                    return;

                inputs = outputs;
            }
        }

        public double testGMDH(double[] x)
        {
            double[] inputs = new double[m];
            for (int i = 0; i < m; i++)
                inputs[i] = x[i];
            int L;
            double[] outputs;

            // формирование рядов
            for (int it = 0; it < rows.Count; it++)
            {
                // формируем вектор выходов
                L = rows[it].L;
                outputs = new double[L];
                for (int l = 0; l < L; l++)
                {
                    double xi = inputs[rows[it].models[l].input1];
                    double xj = inputs[rows[it].models[l].input2];
                    outputs[l] = rows[it].models[l].f(xi, xj);
                }
                inputs = outputs;
            }

            double t_val = inputs[0];
            return t_val;
        }

        void swapModelsInArray(Model[] array, int position1, int position2)
        {
            Model temp = array[position1];
            array[position1] = array[position2];
            array[position2] = temp;
        }
    }
}
